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The past decade has demonstrated increasing interests in using optimal control based methods 
within coherent quantum controllable systems. The versatility of such methods has been demon- 
strated with particular elegance within nuclear magnetic resonance (NMR) where natural separation 
between coherent and dissipative spin dynamics processes has enabled coherent quantum control 
over long periods of time to shape the experiment to almost ideal adoption to the spin system and 
external manipulations. This has led to new design principles as well as powerful new experimental 
methods within magnetic resonance imaging, liquid-state and solid-state NMR spectroscopy. For 
this development to continue and expand, it is crucially important to constantly improve the un- 
derlying numerical algorithms to provide numerical solutions which are optimally compatible with 
implementation on current instrumentation and at same time are numerically stable and offer fast 
monotonic convergence towards the target. Addressing such aims, we here present a smoothing 
monotonically convergent algorithm for pulse sequence design in magnetic resonance which with 
improved optimization stability lead to smooth pulse sequence easier to implement experimentally 
and potentially understand within the analytical framework of modern NMR spectroscopy. 
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I. INTRODUCTION 



Optimal control theory (OCT) is a powerful method for control and design of processes within quantum dynamics. 
Originally the method was applied for problems within engineering and economics. Q, 0| During the past decade or 
so, optimal-control-based methods have been increasingly used for a development of new experiments within optical 



spectroscopy, 
spectroscopy, 



3 



8] 
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quantum information processing, 9Ml4 1 



3 



magnetic resonance imagin g (M RI) 

3 



between electron and nuclear magnetic resonance 
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liquid- and solid-state nuclear magnetic resonance (NMR) 
30-3(1 and dynamic nuclear polarization (DNP) hybrids 



Such applications have not only been useful for the specific 



disciplines taking advantage of new efficient design procedures and improved experimental methods, but it has also 
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47] The latter addresses fundamen- 



stimulated mathematical investigations in quantum optimal control theory, 
tal questions concerned with controllability, convergence, and the establishment of powerful numerical methods for 
optimal control in quantum systems. 

So far, the vast majority of optimal control applications within magnetic resonance have taken advantage of gradient- 
based methods, such as t he g radient ascent pulse engineering (GRAPE) algorithm introduced for NMR applications by 
Khaneja and coworkers, 19( and recently further developed and distributed for general use by Nielsen and coworkers 



in an optimal control version 



21 



29j of the open-source NMR simulation software SIMPSON 
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In combination 



with conjugated gradient algorithms, this method 
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29] proves to be very powerful, as demonstrated by 



numerous applications in which NMR methods with improved experimental sensitivity, robustness towards variations 
in instrumental or spin system parameters, reduction of the undesirable effects from dissipative processes (i.e., re- 
laxation), and experiments with lower radio- frequency (rf) power requirements and thereby reduced risks for sample 



heating have been developed. More recently, it has been demonstrated [38| that a monotonic convergent variant 
of optimal control method based on the algorithms of Krotov 2|-|4j, |6( represents an interesting alternative to the 
gradient-based approaches for efficient experiment design. Our initial work with this algorithm, in a density oper- 
ator formulation, exposed important computational properties of optimal control methods such as global extremum 
searching, fast convergence, algorithmic simplicity, and independence on time discretization in terms of convergence. 

Despite increasing use, it is apparent that current methods face serious challenges in the practical realization which 
needs to be addressed to exploit the full potential of optimal control based methods for design of optimal experiments. 
This applies not least for the most challenging purposes involving large spin systems, optimizations ensuring broadband 
or band-selective operation with respect to certain spin system parameters (e.g., chemical shifts), and powder samples 



in solid-state NMR spectroscopy. For example, looking at the many optimal control pulse sequences proposed so far, 



it ap pea rs that many of these display quite wild oscillations in phase and amplitude of the rf control fields (see, e.g., 
Refs. 



18 



38} These problems have 



2CH29J) which may complicate implementation on available instrumentation with limitations on the speed 
and accuracy of phase and amplitude switching. Furthermore, it turns out that GRAPE displays a quite strong 
dependence on the initial guess of the pulse sequence, dependence on the applied time discretization (i.e., the number 
of pulse variables, and their duration), as well as unpredictable convergence to local extremum points. Along the same 
lines, the monotonic Krotov-based algorithm faces problems with numerical instability and difficulties in a selection 
of parameters controlling the flow of operations and balance of necessary running costs, 
introduced undcsircd needs for intuition, experience, and repetition of optimizations with a very large set of different 
initial guesses in the usage of these methods as replacement for stronger mathematical recipes. Although part of 
these problems have been overpassed in recent adoptions for wave function formalism in optical spectroscopy, |4J| a 
strong need for solutions to the problems still exists for magnetic resonance applications typically performed within 
a density operator formalism. 

In this paper, we present a modified monotonic algorithm which stabilizes convergence and smooths the resulting 
NMR pulses sequences through the use of a frequency truncation technique in course of the optimization. The latter 
aspect is practically important realizing that most optimal control sequences presented so far within NMR spectroscopy 
display quite significant and fast amplitude and phase modulations, which may cause unnecessary problems upon 
implementation on available spectrometer hardware. This work builds on related techniques introduced in the field 
of the laser control of alignment and rotation [50( for an unique controlling field. 

II. THE OPTIMAL CONTROL PROBLEM 

The most typical setup for optimal control pulse sequence design in NMR spectroscopy involves systematic gen- 
eration of optimal radio-frequency (rf) pulse sequences which in a given spin system either (i) accomplish the most 
efficient transfer of coherence or polarization from a given initial spin state po to a desired target spin state C (often 
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52] emphasizing 



referred to as state-to-state transfer) or (ii) synthesize a specific effective (or average) Hamiltonian 
or suppressing certain parts of the internal nuclear spin interactions to tailor the Hamiltonian for evolution under 
desired interactions. These interactions may provide a desired state-to-state coherence/polarization transfer, specific 
information, or provide spectral simplification, e.g., in the form of improved spectral resolution. 

In absence of dissipative processes, the dynamics of the nuclear spin system may be described by the Liouville-von 



Neumann equation 

dp(t) 



= -i[H(t),p(t)} ) (1) 



df 

where p(t) is the density matrix (initial state: p(0) = po) an d H{t) is a Hamiltonian of the spin system. In the 
high-field approximation, the latter takes the form 

H{t)=H +J2"k(t)H k , (2) 

k 

with the first term collecting all internal nuclear spin interactions (chemical shifts as well as J, dipole-dipole, and 
quadrupole couplings) and the latter describing external rf manipulations with the amplitude Wfc(i) (i n angular 
frequency units) for the spin operator H k (typically H k — I x . I y , S x , S y for an I-S two-spin system) being our control 
fields. The solution to the equation of motion in Eq. (TT]) is typically expressed as 

p(t) = U(t)p tf(t), (3) 

where the unitary operator (or propagator) U(t) = D exp (-i f* ff(t')dt') links the unitary evolution with the Hamil- 
tonian in Eq. D is the so-called Dyson time-ordering operator. 
Optimal control relies on optimization of a functional of the type 

Ji(w)=&- VA fc f u a k (t)dt, (4) 
k Jo 

with the first term denoting the final cost (or the objective) and the latter term the running cost considering the 
collected energy /power used to reach the optimum. The influence of the running cost is scaled by a so-called penalty 
factor Afe (which may be constant, as assumed here, or time-dependent if, e.g., specific rise and fall time behaviour of 
the pulses arc desired). This convex running cost improves the convergence of the optimization methods. Furthermore, 
it facilitates development of pulse sequences without too excessive rf power consumption. 

For state-to-state transfers between Hermitian operators C and po, the final cost (i.e., the overall transfer efficiency) 
may be expressed as 

cf> 1 = Tr(Cp(T)), (5) 
where T is the overall duration of the experiment. For transfers between non-Hermitian operators the final cost may 



be formulated as 



29 



2 = |Tr(GV(T))| 2 . (6) 



For synthesis of a desired propagator Ud , the final cost is given by 



3 = Re pR(tf(T)u£) 



(7) 



We note that the objectives <f>i given above is by no means exclusive, and other variants for the target functions may 
be used. For example, instead of Eq.(|7|) one may use the squared expression (f>' 3 = \Tr(u\jJ (T)\ 2 . The preferential 



form of the target function typically depends strongly on the given optimization problem. 
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III. MONOTONIC CONVERGENCE 

The aim of optimal control experiment design is for a given period of time T and a given time discretization (i.e., 
number of pulses, typically taken equidistantly over the time T) to find rf control fields (i.e., ujk(t)) which lead to the 
maximum of the functional in Eq. The Krotov-based monotonic algorithms accomplish this using a Lagrange 

approach with an adjoint propagator B{t), or a Lagrange multiplier for an unconstrained functional as described 



in Ref. ( 38[), using a combination of forward propagation with U(t) (Eq.Q ) and backward propagation using a 
conjugated equation for B{t). We note that different formulations exist for such calculations, including that of Tannor 
and coworkers [3] closely following the original formulations of Krotov, a different variant by Zhu et al. [4], and a 
more general approach (embra cing the methods of Tannor and Zhu) as described by Maday et al. 6] in context of 
waveforms and Maximov et al. |38| in context of density operators. For simplicity, we here restrict to the original 



Krotov approach as formulated by Tannor et al.[3j (corresponding to 5 = 1 and rj = in Refs. ( 6]) and ( 38 
while noting that the more general formulation with arbitrary 5 and rj parameters may readily be implemented at the 
expense of slightly more complicated formulas. 

A prerequisite for monotonic convergence is that the target operator C is positive semi-definite. While this condition 
is obviously fulfilled for optimization of effective Hamiltonians, it will, for example, for optimization of transfer 
between Hermitian operators require the final cost function to be modified to (f> = <ft + kTt ([/(T)[/^(T)) as proposed 



previously. 



38] 



The proof for monotonic convergence, in the case of Hermitian operators, may readily be established considering 
the following decomposition of the variation in the cost functional values between arbitrary rf pulses controls w and 



u/ (in a discrete representation): 

J(w') - J(w) = Tr ((^ - E/*)(E# - + Tr (c(U' N - Efa)po(E# - U N j) 

N+l 

- AtA^^^.! - w k j-x){u/ kij _ x +u> k ,j-i), (8) 

where we have applied the second-order Strang method to evolve the propagator Uj and the Lagrange multiplier 
Bj. 44j The time step At is defined through the number of bins N and the overall time T — NAt, while the matrix 
exponent A is defined as A = e -°- 5lAtH o _ _g j s the identity operator. We assume all functions to be constant during 
the time step At rendering computation of the matrix exponents straightforward. 

The discrete set of propagators (Uj) and Lagrange multipliers ( Bj) may be obtained using the following equations 

Uj+i = Ae- iAt ^" Uk i Hh AUj, 

U = E, (9) 
Bj = B j+1 Ae- iAt ^* Uk i Hh A, 

B N = 0(C), (10) 

with the operator 6(C) defined by Eqs. ©-([7J and the desired efficiency. For example, for Hermitian state-to-state 
transfer 

6(C) = kU n + CU NPo , (11) 



where k is an insignificant scaling factor. 



To ensure positiveness of the functional difference J(oj') — J(w) at each time step, we need to maximize the last string 



4J] it appears convenient to re-express the 



of Eq.®. Following the approach previously described for wave functions, 
two last terms of Eq.® as a function fj(ujj) depending on a vector variable uj'a of dimension k relating to the time 
step j, 

fj(Lj'j) = 2Re ^Tr^(e iAt ^ Uk 'i- lHk e- iAt ^<J-^ H >' - AU'^B^A*} 

- M^Kj^ - ^k,j-i)(^' k ,j-i (12) 
k 

In this formulation, the updated control fields uj' k ■ should be found locally (in time) through a minimization problem 
of dimension k for the vector function —fj(uj^) (corresponding to maximization of fj(ui'j) and the corresponding 
functional in Eq.([5])). which due to the non-commutative relationships of the operators Hk in the matrix exponents 



can not straightforwardly be simplified further. 6 



44j The minimization may be conducted starting with an appropriate 



initial guess, using routines such as conjugated gradients or quasi-Newton. 53] The iteration equation for the rf controls 
uj'j may be expressed as 

wf™ = axgmin {-/><)}, (13) 



with the rf fields L)j from a previous iteration step used as an initial guess for the minimization problem Eq. (|13p . 

Equation (|13|) guarantees monotonic convergence of the functional J (to) for all time steps in each iteration step 
through solution of an unconstrained nonlinear optimization problem. Indeed, fj(^'j = u>j) = 0, so that — /(w™ 6 ™ 1 ) < 0. 
In other words, there exist a pulse sequence uj™ ew such that J (a;™ 6 ™) < J(uij). Moreover, note that if the algorithm 
finds that J(uj™ ew ) = J(ojj) then every local optimization procedure has failed to find a strictly better control. In 
this case, a local maximum has been found. Note also that the previous approach is independent of the optimization 
method used in the local optimization problems in Eq . (fT3l) . 



IV. FREQUENCY CONSTRAINING AND SMOOTHING 



Direct application of the monotonic algorithm outlined above, as well as previous gradient-based algorithms, often 



20, 



23, 



25 



29 



38 



lead to optimal pulse sequences with significant oscillations in the rf field amplitudes and phases 
40j These variations not only hamper analytical understanding of the function of optimal pulse sequences but may also 
complicate practical implementation on available instrumentation. To stimulate generation of smoother solutions, we 
here demonstrate that the optimal control algorithms may be combined with standard frequency truncation techniques 
with a regularization substep that retains monotonic convergence. 

Given cj, an arbitrary control, suppose that a better control ui new has been obtained, i.e. J(w) < J(ui new ). A way 
to define a smoother improving control cj smooth is to consider an interpolation between uj new (i.e., those typically 
displaying significant oscillations) and a regularized version of it (cf. Ref. [50(). In this way define for all k 



smooth 



(14) 



where a is a regularization parameter and ¥{uj) defines the frequency truncation. The frequency truncation may be 



62| or simply by Fourier transforming the 



accomplished in many ways, e.g., using built-in functions in MATLAB 
pulse sequence, removing high-frequency components, and transforming it back to the time domain using an inverse 
Fourier transformation. 



8 

Following the proofs in the previous section, it is obvious that a = ensure monotonic behaviour and thereby 
J(u;™Q° th ) > J(ui). The same is almost true when a — > as well. Accordingly, an iterative sub-algorithm may be 
established which generates a smooth and monotonic solution: 

1. start from a = 1; 

2. test J(uj smooth ) > J(w)? 

3. if not, decrease a, for example, halving the value and go to step 2. 

With this ingredient and the formula in Eqs. ([9]) - ([TT|) . we can formulate the overall algorithm. Consider a 
parameter Tol > 0. 

1. Set the initial random guess lj° = (cjj k ); j = 1, .., N; k > 0. 

2. Compute the corresponding state Uo and B according to Eqs. (P)- (fTU|) . 

3. Set e = +oo, I = 0. 

4. While e > Tol do 

a. Do forward propagation and search new rf field ui l+1 according to the procedure in Sec. Ill 

b. Apply the frequency truncation sub-algorithm with regularization to find a 1 that preserves the monotonicity, 
i.e.: J((l - a e )u e+1 + a e F{u e+1 )) > J{u: e ) 

c. Define u t+1 = (1 - a l )u l+1 + a l F{u l+1 ) 

d. Sete= J{oj e+1 )- J(c/) 

e. Do backward propagation to compute B l+1 the solution of Eq. (JTUJ) with lu = us e+1 

f. Do i = I + 1 
End While. 

It is important to note that frequency truncation technique with regularization could equally well be applied to 
other known optimal control approaches, such as GRAPE and different variants of Krotov-based optimal control. 

V. DESIGN OF SMOOTH NMR EXPERIMENTS 



While the overall applicability of optimal control for NMR experiment design has been demonstrated and verified 
numerically and experimentally in numerous previous papers, 15l - ll7lll9l 422. 24, 26-l29j we will demonstrate numerically 



9 

in this paper the implication of our smoothing procedure and monotonic convergence in optimal control design of a set 
of simple NMR experiments. We address specifically spin-state selective, non-Hermitian coherence transfer through 
J couplings for liquid-state NMR applications and excitation of the central transition of 23 Na for applications in 
magnetic resonance imaging (MRI). Lists containing rf amplitudes and phases for the optimal control pulse sequences 
as well as Matlab scripts used to generate these can be found in Supplementary Material. [63] 



A. Coherence-order and spin-state-selective coherence transfer 

For a two-spin, J-coupled spin systems in the context of liquid-state NMR, the internal Hamiltonian may in the 
high-field approximation be cast as 

H = nJ2I z S z , (15) 

where I z or S z represent z-components of the spin operators for the I and S spins in the present example coupled 
through a scalar coupling of size J = 140 Hz. To exemplify transfer between non-Hermitian operators, we assume 
the initial state to be represented by +l-quantum coherence on the S spin (i.e., po = S + ) and to further demonstrate 



spin-state-selective transfer, 



•54 



- 57] we assume that the destination operator is -1-quantum coherence on the I spin 
with the S spin being in the a-state corresponding to only one of the lines in the J-coupled doublet excite (i.e., 

c=i-s a ). 

In this case of transfer between non-Hermitian operators, we use the target function in Eq.® modified to ensure 
our target being positive semi-definite. This leads to the modified objective <j) and the operator 6(C): 

t = <f>2 + Tr(U(T)tf(T)), (16) 
9(C) = U\T) + Po UHT)C* l Tr \u 1 \T)p\u{T)C 

+ CtC/VoTr [U{T)p U\T)&] , (17) 

where we, relative to the expression in Eq. (jlip . assumed the scaling factor to be k e 1 for the sake of simplicity. 
Our control fields are represented by the amplitudes (angular frequencies) Wfc(i) of x- and y- phase rf irradiation on 
the spins I (operators I x (fc=l) and I y (k=2)) and S (operators S x (k = 3) and S y (k — 4)). According to unitary 



bounds on spin dynamics. 



58, 



59J the maximum achievable transfer efficiency in this case is 1. For the optimizations 
we set the excitation period to T — 7.14 ms (corresponding to 1/J) and the number of pulses to N — 200. 

Addressing this specific optimization problem, Figure Q] compares pulse sequences (i.e., I and S spin rf field 
strengths) obtained on basis of a random initial pulse sequence using the gradient-based optimal control algorithm 



GRAPE. 
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38] and the latter combined 



19j the monotonic Krotov-based optimal control algorithm of Maximov et al., 
with smoothing as described in this paper. It is quite evident that the pulse sequences developed using the two former 
methods display quite significant oscillations, for which the appearance depends very much on the specific optimiza- 
tion. In the present case it looks like the GRAPE algorithm produces a more smooth, albeit still oscillating a lot in 
the "waves", pulse sequence than the Krotov-based algorithm. Slight variation in parameters, such as the number 
of steps in the waves, may radically change this pictures and in a later example we will see the opposite picture. 
The same variability will be seen upon inspection of the many optimal control NMR pulse sequences presented in 
literature so far. What is clear, however, is that the Krotov-based algorithm (as selected in this paper, but it could 
just as well be GRAPE) combined with frequency-truncation smoothing leads to much smoother sequences which may 
easier be analysed analytically and which put less demands on the spectrometer hardware upon practical implemen- 
tation. Apart from this the sequences display only relatively modest variations with respect to rf power consumption 
(root-mean-square (RMS) average rf powers (I and S spin values separated by /) of 141/71 Hz for GRAPE, 92/91 
Hz for Krotov, and 12/102 Hz for smoothing Krotov) and are essentially identical with respect to coherence transfer 
efficiency (99% of the theoretical maximum for all methods). 

While rigid statements on the optimization speed and convergence require a much more thorough analysis (to 
be presented elsewhere), examples like the one in Fig. [T] provide the following crude estimate. Smooth Krotov is 
computation-speed wise quite similar to GRAPE, while the original Krotov approach is somewhat faster (in the order 
of a factor 3-5, although this number highly depends on the specific optimization). Out of 1000 optimizations based 
on random initial pulse sequences only 80-90 % of the GRAPE and Krotov optimizations led to sequences with more 
than 90% of the nominal transfer efficiency, while all sequences in this specific case passed this limit for the smooth 
Krotov approach. The great benefit of smooth Krotov is that it provides smooth pulse sequences, and the benefit of 
both Krotov variants is that they enable optimization with much more coarse time dicretization (longer and fewer 
pulses) than the GRAPE algorithm due to its fundamentally different optimization strategy. 

Figure [2] gives snapshots of pulse sequences throughout optimization using the smoothing Krotov-based optimal 
control algorithm starting out from a random pulse sequence (with maximum amplitude of 100 Hz). The snapshots 
illustrate gradual adoption of the optimal sequence to a smooth appearance. It is evident that already after 5 
iterations, the control fields on the / spin almost vanishes, leaving the major action to the S spin control fields for the 
remaining optimization into the final pulse sequence. The final pulse sequence have quite low rf-power consumption, 
due to the reducing effect of the running cost (Ai = 10 -4 s 2 rad~ 2 , A2 = Ai for both the / and S spins), and the rf 
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fields vary smoothly as an effect of the smoothing algorithm. The progress of the cost function and its constituents 
(the overall functional (cost), the transfer efficiency, and the running cost (penalty)) throughout the optimization is 
illustrated in Fig. [3] which also provide a numerical demonstration of monotonic convergence. It is evident that the 
functional converges to a value which, when considering the subtractive term from the running cost, is equal to or 
close to the theoretical limit after a number of iterations where the efficiency (final cost) and penalty (running cost) 
have displayed some exchanging oscillations before converge to the optimal values. 



B. Optimal control for satellite and central transition excitation in 23 Na MRI 

Optimal control mediated experiment design is by no means restricted to spin-1/2 nuclei or systems of these. 
Many challenging optimization examples may be found for NMR spectroscopy or MRI in concern of quadrupolar 
nuclei, as rece ntly demonstrated by optimal control pulse sequences for improved multiple-quantum magic-angle- 
spinning NMR[60| and for selective excitation of the central transition of 23 Na (spin I = 3/2) in presence of residual 
quadrupole couplings for magnetic resonance imaging purposes. 



34 



The latter application may be very important 



for 3 Na MRI of, e.g., cartilage where the sodium concentration may be considered a reporter for disorders and 



degradation. 



34, 



6l| In such applications, the ability to separate 23 Na ions with large and very small (vanishing) 
residual quadrupolar couplings is regarded important as they represent different populations of relevant ions. 

In this section, we demonstrate the use of the smoothing Krotov-based optimal control algorithm to design pulse 
sequences which selective excites the central transition (— or the satellite (— §,— i or i,|) transitions of 23 Na ions 
characterized by a residual quadrupole coupling frequency of wq/27t = 60 Hz. In this case the size of the quadrupole 
coupling and the rf irradiation fields are comparable (often referred to as the intermediate regime) implying that 
experiment design by standard analytical means is not straightforward. 

The optimization involves the secular, high-field-approximated first-order Hamiltonian 

H=^-(3^-I{I + l))+u 1 {t)I x +u 3 (t)I y , (18) 

where ujq is quadrupole coupling frequency while LO\{i) and uJ2(t) are control rf fields corresponding to the spin 
operators I x and I y , respectively (all frequencies in angular units). The nature of the optimization will be a standard 
state-to-state transfer with the final cost expressed by Eq. ([5]) with the destination operator C representing x-phase 
coherence on either the central transition or the satellite transitions. The initial state corresponds to the thermal 
equilibrium p = I z . 
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Figure [4] illustrates the performance of optimal control pulse sequences designed to enhance the central transition 
(in this case leading to efficient suppression of the satellite transitions although not requested in the cost function 
used for the optimization) or excite the satellite transitions while suppressing the central transition (also here the 
suppression was not specified in the cost function) relative to a hard (infinitely strong) 90° pulse exciting both types 
of transitions. The simulated spectra clearly reveal that (i) the central-transition selective excitation optimal control 
sequence offers a sensitivity enhancement, relative to standard single-pulse excitation, by the maximal factor of 1.5 as 



,3 



discussed previously by Jerschow and coworkers, |32j while ensuring reasonable suppression of the satellite transitions 
and (ii) the pulse sequence designed to suppress the central transition lead to satellites with unaltered intensity while 
the central transition is almost removed. 

The optimal control pulse sequences providing selective excitation of the central transition and the satellite tran- 
sitions are shown in Figs. [SJ\ and[5|3, respectively. The sequences have a overall length of 2.25//q (with the peak 
splitting fq related to the quadrupole coupling as fq — Zujq/2n) and accommodates 200 pulses of equal length. The 
central transition excitation pulse sequence is characterized by a RMS average rf power of 45 Hz and excites the cen- 
tral transition with an efficiency of 99% and the satellites have been suppressed to an efficiency of 1.4% (percentages 
relative to the theoretical maximum). The satellite transition sequence is characterized by an RMS rf power of 54 
Hz, with the satellites excited with an efficiency of 99% and the central transition suppressed 72%. 

It is clear from both sequences that the optimal pulse sequence is found in the regime u> r f ~ uiq for which experiment 
design by standard analytical means is exceedingly difficult. It is also evident that the smoothing algorithm used 
during optimization leads to smooth sequences offering easier implementation on conventional MRI instrumentation, 
than the corresponding sequences we obtained using previous optimal control formulations. The latter aspect is 
demonstrated in Fig. [BJ'V by central transition-selective excitation sequences (corresponding to Fig. obtained 
using GRAPE or in Fig. [6j3 using our previous Krotov-based optimal control software. The latter pulse sequences 
involved 500 pulses for GRAPE and 200 pulses for the Krotov-based approach to facilitate comparison with the 
pulse sequences earlier proposed by Jerschow and coworkers [34j for GRAPE. In the present implementation, the 
GRAPE-derived sequence displays much wilder oscillations than the Krotov-based sequence. 



C. Broadband optimization 



In practical applications of optimal control theory for experiment design, it is often desirable to optimize the pulse 
sequences to be robust toward variation in one or more spin system parameters, introducing the concepts of broadband 
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24j and Krotov-based 38[ methods 



or bandselective excitation as described previously in the context of GRAPE 
Such requests may readily be incorporated into our optimal control algorithm through a modified cost function 

J(w) = kTy(U(T)U\T)) +J2 r ^(C 1 U i (T)p UJ(T)) 

i 

- \J2 f "l(t)dt, (19) 

where each propagator Ui{T) corresponds to a specific condition, such as different chemical shifts, couplings, or 
orientations of the sample relative to the magnetic field (in solid-state NMR), where U is a common propagator of 
the system. The additional targets in the functional Eq. (|19[) corresponding to the broadband extension demands 
a modification of the condition of a monotonic convergence and as concequence changing the form of Eq. (1121) . This 
leads to a modification of the target function fj(u>'j) in Eq. (ll2l) where the part involving propagators and Lagrange 
matrices should be replaced by a summation 

••• Y, AiU'^B^Al... (20) 
conditions i 

with the propagator Ui and adjoint matrix Bi computed for each individual condition when the pulse sequence ujk is 
the same for all propagators. 

To demonstrate the application of broadband optimization, we repeat our optimization of 23 Na MRI central tran- 
sition enhancement experiments for a range of quadrupole coupling around ujq/2tt = 60 Hz with a deviation of ±20 
Hz. This leads to a pulse sequence with a central transition excitation profile as function of the quadrupole coupling 
frequency as illustrated in Fig. [7J revealing very good excitation in proximity of ujq/2tt — 60 Hz as requested in 
the optimization. We note that it is obviously possible to specify specific regions of excitation and no excitations by 
entering additional conditions in the target function and algorithm. 



VI. CONCLUSION 



In conclusion, we have presented a smoothing monotonically convergent optimal control algorithm for efficient 
design of pulse sequences for magnetic resonance applications. The proposed frequency-truncation algorithm has 
been incorporated into a Krotov-based optimal control procedure and demonstrated for a few of NMR examples 
where it leads to much smoother optimal pulse sequences than obtained using previous methods. We note that 
previous methods may eventually lead to smooth sequence by themselves, in particular for simple optimizations 
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without inhomogeneities and broadband performance request and in cases where educated guesses are present to 
initiate the optimization. However, in the most general case optimal control typically leads to pulse sequences with 
significant oscillations in phase and amplitude. In such cases our new algorithm offers a remedy to incorporate 
smoothness in the optimization. The smooth sequences improves the possibilities for obtaining theoretical insight 
for the optimal sequences and facilitates implementation on typical spectrometer hardware. We anticipate that the 
methods will find widespread applications for design of experiments within liquid- and solid-state NMR spectroscopy, 
MM, and dynamic nuclear polarization. 
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VII. FIGURE CAPTIONS 



Fig. 1. (Color online) Pulse sequences for spin-state- and coherence-order-selective S + — > I~S a coherence transfer 
in a J-coupled (J=140 Hz) heteronuclear two-spin system designed using optimal control based on the GRAPE 
algorithm (first row), the original monotonic Krotov-based algorithm (second row), and the latter combined with 
frequency-truncation smoothing (third row). The red and green lines correspond to x- and y-phase rf control fields 
for the spins I (left column) and S (right column), respectively. 

Fig. 2. (Color online) Representative pulse sequences obtained throughout the optimization in Fig. 1, illustrating 
the gradual adoption of the pulse sequence to provide optimum S + — > I~S a coherence transfer with modest rf-power 
consumption and smooth rf variation. The red and green lines correspond to x- and y-phase control fields for spin I 
while blue and pink lines represent x- and j/-phases of the 5-spin control fields. 

Fig. 3. (Color online) The dependence of the functional J{uj), the efficiency <j) 2 , and the running cost (penalty) 
throughout the iteration steps of a smoothing monotonic Krotov-based optimal control design of the pulse sequence 
for S + — )■ I~S a coherence transfer, shown and analysed in Figs. 1 and 2. The red curve represents the functional 
J(co) (defined by Eqs. (4) and (16); for transparency, we have subtracted the correction terms in Eq. (16) which 
otherwise would add a constant contribution to the functional in all points), the green the efficiency <f)2, and blue the 
penalty cost (see text). All ordinates have been scaled by a factor Tr{CC^}. 

Fig. 4. (Color online) Simulated 23 Na NMR spectra of sodium with a residual quadrupole coupling frequency of 
loq/2tt — 60 Hz excited using a standard 90° hard pulse (red line), the optimal control pulse sequence in Fig. 5 A 
designed to excite exclusively the central transition (green line), and the pulse sequence in Fig. 5B designed to excite 
the satellite transitions only (blue line). 

Fig. 5. (Color online) Optimal control pulse sequences for 23 Na MRI, which for the case of a residual quadrupole 
coupling of 60 Hz offer selective and enhanced excitation of the central transition (A) or excitation of the satellite 
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transition while suppressing the central peak of the 23 Na NMR/MRI spectrum (B). The red curve represents x-phase 
control fields while the green curve represents y-phase controls. The length of the pulse sequences are 12.5 ms. 

Fig. 6. (Color online) Optimal control pulse sequence for excitation of the central peak of the quadrupole spectrum 
of 23 Na (A) (parameters as described in relation to Fig. 4) obtained using the GRAPE algorithm and (B) the 
original Krotov based algorithm. The red and green curves represent x- and t/-phase control fields. The length of 
pulse sequences in both cases is 12.5ms. 

Fig. 7. (Color online) (a) Excitation profile (normalized) for the central transition as function of the quadrupole 
coupling frequency for a pulse sequence optimized using our smooth variant of the Krotov algorithm and specifying 
"broadband" excitation for a range of quadrupole coupling frequencies in proximity of 60 Hz. (b) The length of the 
pulse sequence is 12.5 ms, the RMS average rf power is 170 Hz, and the average transfer efficiency in the range of 
quadrupolar couplings between 40 and 80 Hz is 1.89 (corresponding to 95 % of the theoretical maximum). 
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